Entrainment transition in populations of random frequency oscillators 
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The entrainment transition of coupled random frequency oscillators is revisited. The Kuramoto 
model (global coupling) is shown to exhibit unusual sample-dependent finite size effects leading to a 
correlation size exponent v = 5/2. Simulations of locally-coupled oscillators in d dimensions reveal 
two types of frequency entrainment: mean- field behavior at d > 4, and aggregation of compact 
synchronized domains in three and four dimensions. In the latter case, scaling arguments yield a 
correlation length exponent v = 2/(d — 2), in good agreement with numerical results. 
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Collective oscillations are abundant in physical, chemi- 
cal, and biological systems far from equilibrium 

[EH Si, 

5] . Such synchronized behavior has been widely explored 
via various systems of coupled oscillators [a 0, B E3] • 
Although several theoretical methods are available to 
treat systems of identical oscillators with or without 



noise 

lid, 

I 111 . |X2I ] . the description of synchronization in 
coupled oscillators with a broad distribution of intrinsic 
frequencies is still a largely unsolved problem [g, |9( . More 
recently, the study of dynamics on complex networks, in- 
cluding the human brain, brought renewed interest into 
the subject [ijj]. 

The onset of system- wide synchronization at some cou- 
pling strength is often compared to critical phenomena 
in equilibrium systems, where scaling concepts provide 
an adequate theoretical basis for quantitative analysis. 
The scaling approach has been successfully applied to 
the Kuramoto model but its extension to finite- 
dimcnsional systems has stumbled over a number of diffi- 
culties 0, [II [3 El. Even in the globally-coupled case, 
the critical size needed to maintain a coherent group of 
synchronized oscillators in a large but finite population, 
as well as their ensuing dynamics, have not been fully un- 
derstood d, @, [TtJ ■ In d-dimensional systems with local 
coupling, d = 2 is generally accepted as the lower critical 
dimension for macroscopic entrainment, but the upper 
critical dimension is uncertain due to peculiar size effects 



17] 



seen in numerical investigations 

In this Letter, we re-examine the above issues criti- 
cally and establish two types of behavior at the onset 
of macroscopic synchronization in coupled random fre- 
quency oscillators. In the globally coupled case and in 
finite dimensions d > 4, we show that frequency entrain- 
ment is accompanied by spontaneous symmetry break- 
ing in the phase, and as such shares many features with 
the usual ordering transition in equilibrium systems with 
0(2) symmetry. The correlation size exponent v gov- 
erning finite-size scaling at criticality, however, takes on 
the value i in the present case rather than 2 in usual 



mean- field theories [18|. For d < 4, the above type of 
symmetry-breaking is ruled out by diverging phase fluc- 
tuations. Instead, oscillators are first entrained with their 
neighbors to form compact, synchronized domains which 
grow upon further increase in coupling. System-wide 
entrainment is achieved only when the largest domain 
reaches system size. This gives rise to size- and sample- 
dependent thresholds whose behavior is quantified. 

We begin with the globally-coupled Kuramoto model 
defined by the following set of dynamical equations: 



d<f)j 
~dt 



KAsmfa - 9). 



Here <f>j is the phase of the jth oscillator (j = 1, 



(1) 

,N), 

and ojj its natural frequency, drawn from a normalized 
distribution g(u>). To be definite, here we take g(of) to 
be Gaussian with zero mean and unit variance. The pa- 
rameter K sets the strength of coupling to the global 
quantities A and 9 defined via 



Ae 



N' 



N 

J'=l 



(2) 



In the classic work by Kuramoto, the existence of an 
entrained state is established by considering solutions to 
(H|) at a constant A. After an initial transient, oscillators 
with \u)j | < KA reach a fixed angle while those with 
\u>j\ > KA are in a running pendulum state. In the limit 
N — ► oo, the self-consistent equation for A reads 



KA 



dug(u))Jl 



-KA 



( — ) 2 



*(A) 



(3) 



which has a nontrivial solution when K > K c 
For a population of finite size, 'J (A) is replaced by 



3,\uj\<KA 



l-(i^) 



7rg(0) ' 



(4) 



2 



which contains a sample-dependent correction 5^> = — 
(vj>) = f - $ oc y/W^/N. Here TV S is the number of 
oscillators in the frequency interval (— K A, ATA), and (•) 
denotes sample average. Close to the transition, the self- 
consistency equation A = ^(A) can be written as 



A = aKA - c(KAf 



5\T>, 



(5) 



where a = K c 1 and c = — 7rg"(0)/16. The variance of 8^> 
is given by ((S^) 2 } = \g{Q)KA/N + 0{A 2 /N). Hence 
the solution to Eq. §5§ takes the scaling form 



A(K, TV) = N- 1/5 f(kN 2/b ), 



f 2/ 5) 



(6) 



where k = (K — K c )/K c . The scaling function f(x) is 
sample-dependent and satisfies the equation, 



xf-cK 3 c f + (8/3^ 2 ef^ = 0. 



(7) 



Here e = 5S& / {(5S&) 2 ) 1 ' 2 is a Gaussian random variable 
with zero mean and unit variance. 

Equation (|6|) resembles the finite-size scaling form in 
various mean-field models of equilibrium phase transi- 
tions [lij], but the value of the exponent D = | is quite 
unusual. To check that A in the transition region is 
dominated by S^> arising from "density fluctuations" of 
oscillators along the frequency axis, we have performed 
extensive simulations of Eq. (TTJ). Figure QTa) illustrates 
sample-to-sample fluctuations in the time-averaged value 
of A 2 . The onset of entrainment spreads over a dis- 
tance 8 K ~ TV~ 2 / 5 around the nominal K c , as seen from 
the scaling plot of (A 2 ) against K at various sizes N in 
Fig. [Jib) . The dashed line there is obtained by averaging 
solutions to Eq. ([7]) at many different values of e drawn 
from a Gaussian distribution. The agreement between 
the predicted scaling and simulation data is satisfactory 
on the entrained side and sufficiently close to the transi- 
tion, on which Eq. ((5|) is based. FigureQJc) shows the rel- 
ative strength of temporal fluctuations SA(t) = A(t) — A 
of the order parameter in the transition region. (Here 
the ovcrlinc bar denotes time average.) On the detrained 

2 

side, the ratio (<5A 2 )/(A ) at large N approaches the 
value — — 1 for Gaussian fluctuations. On the entrained 
side, this quantity decreases as A -1 as expected from 
independent fluctuations. At K = K c , the data clearly 
decrease with TV, lending further support to our analysis. 

We now turn to locally-coupled oscillators on d- 
dimensional hypercubic lattices, for which the dynamical 
equations take the form 



dt 



E 



,0- 



(8) 



where Aj is the set of all nearest neighbors of site i. 

At sufficiently large K, Ecl JHJ) can be treated using 
a linear approximation [14, Il7|. In this limit, the system 
enters a "fully entrained" state where the random term 




FIG. 1: (Color online) A 2 against K for globally-coupled 
oscillators. (a) Time-averaged value of A 2 for 20 differ- 
ent samples (TV = 12800). (b) Scaling plot for TV = 
100, 200, . . . , 12800, in increasing order as indicted by the ar- 
row (averaged over 100 samples for each TV). Dashed line: 
see text, (c) Relative strength of temporal versus sample-to- 
sample fluctuations of A for a set of system sizes as in (b). 



oji is balanced by local gradients of a static phase field 
. Salient features of this state are summarized as 
follows: i) For d > 4, <j)f^ has bounded variations even 
when the linear system size L — > oo. Consequently, the 
entrained state has a global phase that breaks the 0(2) 
symmetry, ii) For 2 < d < 4, variations in 4>f^ grow as 
£( 4 - rf )/ 2 ^ Hence an entrained state cannot be assigned a 
definite phase, iii) For d < 2, local phase gradients have 
an infra-red divergence _L( 2-d )/ 2 . In this case, the system 
is detrained beyond a coherence length £ ~ x 2 /( 2 - d ) . 

As K decreases, the local phase gradients need to in- 
crease to counter the w,-'s but there is a limit to how far 
this can go in the Kuramoto model. Upon proliferation of 
phase slips and runaway oscillators, two scenarios can be 
contemplated for the destruction of global entrainment: 

i) Oscillators break away individually from the entrained 
group until the latter is reduced to an infinitesimal frac- 
tion of the system as in the globally-coupled model, ii) 
A more collective form of phase slips takes place along 
"domain boundaries" that break the system into locally 
synchronized clusters, starting from the largest scale. For 

ii) to preempt i), there must be pre-existing large phase 
differences (i.e., "strain") across the system, which is the 
case in the entrained state at and below d = 4 but not 
above. This suggests that the nature of the detrainment 
transition is quite different above and below d — 4. 

To characterize system- wide coherent phase motion in 
the entrained state, we introduce the Edwards- Anderson 
order parameter 



lim 

-t — K 



1 

TV 



yyw>i(*)- 



</>i(to)] 



(9) 
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Compared to A defined by Eq. |2]| , Aea is incensitive to 
the static phase deformation discussed above, and hence 
can also be used for d < 4. 

We have carried out extensive numerical simulations 
of (H|) to investigate the two types of transition behavior 
and associated critical properties and finite-size effects. 
Periodic boundary conditions are used. Let us first ex- 
amine the spatio-temporal behavior of the phase advance 
A<pi = <f>i(t + T) — 4>i (t) of oscillators over a sufficiently 
long interval T in a given oscillator population near the 
entrainment transition. Figure HJa) presents the distri- 
bution P(v) of the mean phase velocity Vi = Afa/T over 
an interval T = 10000 for a JV = 16 5 system in d = 5 
dimensions. The inset shows a magnified plot of the peak 
region where values of A<j) within each 2ir interval are re- 
solved. Evidently, entrainment here is accompanied by 
the selection of a global phase, i.e., the transition is of 
the symmetry-breaking type. Comparing the two distri- 
butions at K = 0.2 and K — 0.205, we see that only 
a small fraction of oscillators participate in the onset of 
entrainment, while the rest are not much affected at this 
stage. This is very similar to the behavior of the globally- 
coupled case. At K = 0.205, the wings can be fitted well 
to a weak, integrable power-law, P(v) ~ |w| -1 / 2 . 

Figure EJb) shows P(v) for a d = 3 system of N = 64 3 
oscillators at three different values of K around the en- 
trainment transition. Here T = 5000. A strong narrow- 
ing of P(v) is seen in the entire critical region, indicating 
the formation of large synchronized domains well below 
global entrainment. The wings fall off roughly as v~ 2 . 
Consequently, most oscillators are moving at phase ve- 
locities close to that of the peak. However, as seen in the 
inset, the actual phase advance Acpi of these oscillators 
is much less entrained as compared to the d = 5 case. 

Further indication of entrainment through growth and 
aggregation of locally synchronized domains is found in 
the spatial structures of the A^'s, as depicted in Fig. [3] 
for one layer of the d = 3 system. At K = 0.66, a 
large, spanning domain of entrained oscillators coexist 
with smaller clusters of oscillators at varying phase ve- 
locities. When K is decreased to 0.64, part of this largest 
domain undergoes a phase slip over the time period, as in- 
dicated by the arrow in Fig.[3^b). This process of detrain- 
ment through phase slips continues down to the smallest 
scale upon further weakening of the coupling. 

The actual dynamics of phase slip initiation and prop- 
agation, particularly in the presence of runaway oscil- 
lators, is rather complicated and will be left for future 
investigation. It is, however, possible to gain some in- 
tuition about the factors governing the typical domain 
size through the following consideration. In Fig. [3jb), 
for example, each of the light-colored regions can be con- 
sidered as critical in the sense that a weaker K would 
lead to fragmentation of the domain while a stronger 
K would lead to merging with its neighbors. Let SK 
be the increment in K needed for the latter process to 
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FIG. 2: (Color online) Distribution of phase velocities at var- 
ious values of the coupling strength K in (a) five, (b) three 
dimensions. Insets show the peak region in detail. The time- 
averaged value of the order parameter are listed. 



(a) K = 0.62 (b) K = 0.64 (c) K = 0.66 




FIG. 3: (Color online) Spatial structure of <f>i(t + T) - 4> t {t) 
in one layer from the sample depicted in Fig. [2Jb). Here 
T = 5000. Color scale from blue: -12tt to red: 12tt. 



occur when the typical size of light-colored domains is 
£. Upon synchronization of two neighboring domains, 
their phase difference A0o ~ ^( 4 ~ d )/ 2 needs to be accom- 
modated. This is particularly so at the slip boundary, 
where bonds are turned from being barely unstable to 
barely stable. The strength of these bonds is of order 
SK, which offsets the phase gradient A</> /£ ~ £( 2 ~ d )/ 2 . 
Let K c be the value of K when £ = oo, the above analysis 
yields a prediction for the size of synchronized domains 
\ ~ 5K- V ~ (K c - K)' y at a given K, where v = 

Figure [4] shows the entrainment order parameter 
against K for various system sizes. Thirty to several 
hundred samples were used to obtain the average in each 
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FIG. 4: (Color online) Order parameter against K for various 
sizes (L values given on each graph) in three to six dimensions. 
Values of K c [indicated by arrows in (a)-(d)] and v used in the 
scaling plots are (e) K c = 0.8, u = 2; (f) K c = 0.318, f = 1; 
(g) K c = 0.203, v=\\ (h) K c = 0.158, v = 

case. To verify the scaling predictions given above, we 
have analyzed the data with the help of the usual finite- 
size scaling ansatz, 

(A 2 ) = L- 2f3/v <S>{kL 1/u ), (10) 

where k — [K — K c )/K c and [3 and v are the order pa- 
rameter and correlation length exponents, respectively. 
For d = 5 and 6, excellent data collapse is achieved using 
/? = | and v = K = ^ as in the globally-coupled case. 
Interestingly, the scaling extends well into the detrained 
phase. Noting that \ — N(A 2 ) corresponds to the sus- 
ceptibility in equilibrium magnetic systems, we conclude 
that the exponent 7 = dv — 2/3 = | for d > 4 and differs 
from the globally coupled case fl9l | . 

For d = 3 and 4, there is no convergence with increas- 
ing system size (left panel) until (A| A ) reaches a value 
near one. This behavior is consistent with the idea that 
global entrainment takes place only when locally synchro- 
nized domains join to span the whole system. Overall, 
a good data collapse is achieved in each case using the 
predicted exponents = and v = _ 

In summary, through analytical arguments and large 
scale simulations of the Kuramoto model in finite dimen- 
sions, we have established two types of critical behavior 
at the onset of global entrainment. Above four dimen- 
sions, entrainment breaks the global phase symmetry, as 



in the globally-coupled model, with identical scaling ex- 
ponents /3 = I and v = |. For 2 < d < 4, and in par- 
ticular for the physical dimension d = 3, global entrain- 
ment (detrainment) occurs via the aggregation (fragmen- 
tation) of synchronized domains. The size of such do- 
mains obeys scaling with an exponent v = 335. 
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